
Module QMedian

! Quick median by Wirth: Algorith + data structures = programs

  implicit none

contains

  subroutine QMed(n,a,k,x)

    integer, intent(in) :: n,k
    real :: a(:),w,x
    integer :: l,r,i,j

    l = 1
    r = n
    do while( l < r )
       x = a(k)
       i = l
       j = r
       do
          do while( a(i) < x )
             i = i + 1
          enddo
          do while( x < a(j) )
             j = j - 1
          enddo
          if( i <= j ) then
             w = a(i)
             a(i) = a(j)
             a(j) = w
             i = i + 1
             j = j - 1
          endif
          if( i > j ) exit
       enddo
       if( j < k ) l = i
       if( k < i ) r = j
    enddo

end subroutine

function select(arr,k)

  ! select by numerical recipes in F90

  integer, intent(in) :: k
  real, dimension(:), intent(inout) :: arr
  real :: select

  integer :: i,r,j,l,n
  real :: a, b

  n = size(arr)
  if( .not. (1 <= k .and. k <= n) ) then
     select = huge(select)
     return
  endif
  l = 1
  r = n
  do
     if( r-l <= 1 ) then
        if( r-l == 1 .and. arr(l) > arr(r) ) then
           a = arr(l)
           arr(l) = arr(r)
           arr(r) = a
        end if
        select = arr(k)
        return
     else
        i = (l + r)/2

        ! arr(i) <-> arr(l+1)
        a = arr(i)
        arr(i) = arr(l+1)
        arr(l+1) = a
        
        ! arr(l) <-> arr(r), if  arr(l) > arr(r)
        if( arr(l) > arr(r) ) then
           a = arr(l)
           arr(l) = arr(r)
           arr(r) = a
        endif
        
        ! arr(l+1) <-> arr(r), if  arr(l+1) > arr(r)
        if( arr(l+1) > arr(r) ) then
           a = arr(l+1)
           arr(l+1) = arr(r)
           arr(r) = a
        endif
  
        ! arr(l) <-> arr(l+1), if  arr(l) > arr(l+1)
        if( arr(l) > arr(l+1) ) then
           a = arr(l)
           arr(l) = arr(l+1)
           arr(l+1) = a
        endif

        i = l + 1
        j = r
        a = arr(l+1)
        do
           do
              i = i + 1
              if( arr(i) >= a ) exit
           enddo
           do
              j = j - 1
              if( arr(j) <= a ) exit
           end do
           if( j < i ) exit
           
           ! arr(i) <-> arr(j)
           b = arr(i)
           arr(i) = arr(j)
           arr(j) = b

        end do

        arr(l+1) = arr(j)
        arr(j) = a
        if( j >= k ) r = j - 1
        if( j <= k ) l = i
     endif
  enddo

end function select

end module



